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Theory of Impurity Effects on the Spin Nematic State 
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The effect of magnetic bond disorder in otherwise antiferro nematic ordered system is in- 
vestigated. We introduced triangular-shaped ferromagnetic bond disorder in the S=l bilinear- 
biquadratic model on a triangular lattice. It is shown that the coupling between the impurity 
magnetic moment and nonmagnetic excitation in the bulk yields single-moment anisotropy and 
long-range anisotropic interaction between impurity magnetic moments. This interaction can 
induce unconventional spin- freezing phenomena observed in triangular magnet, NiGa2S4. 
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1. Introduction 

Diverse novel low-energy behaviors of geometrically 
frustrated magnets have attracted much attention. ^ 
The central issue is the possibility of a spin liquid, 
namely a quantum disordered state where magnetic 
long-range order is destroyed by frustration and quan- 
tum fluctuation. 2 ) This idea was first proposed by An- 
derson for a Heisenberg antiferromagnet on the trian- 
gular lattice. 3 ) Although subsequent numerical works 
showed the presence of magnetic long range order for 
that model, 4 ) the possibility of spin liquid states have 
been intensively studied, both experimentally and the- 
oretically. Recently various compounds such as organic 
k-(BEDT-TTF) 2 Cu 2 (CN) 3 5 ) and NiGa 2 S 4 6) as well as 
SrCg_ p Ga 1 2-9pOi9 7 - ) are found to exhibit spin-liquid-likc 
behaviors. In addition to their "spin liquid" states, spin 
glass states have been widely observed in geometrically 
frustrated magnets, 1 - 1 and this may be induced by a small 
amount of quenched disorder. These spin glass states 
could be ascribed to coexistence of intrinsic geometrical 
frustration and extrinsic frustration induced by disorder, 
but are not well understood theoretically. Therefore this 
is a challenging problem. Further it is interesting that 
these spin glass states are often observed to accompany 
spin liquid behavior. This implies that these two effects 
are closely related. One typical example of this coexis- 
tence is the case of the triangular lattice antiferromagnet 
NiGa2S4 and this is the issue of this paper. 

In the layered chalcogenide NiGa 2 S4, magnetic Ni ions 
form a perfect triangular lattice well separated by GaS 
polyhedra, and therefore the Ni magnetism has a quasi- 
two-dimensional nature. 6 - 1 Each Ni ion has the electronic 
configuration of t^e 2 and has spin 5 = 1 formally 
with no anisotropy, which is consistent with the nearly 
isotropic susceptibility. Several low-temperature proper- 
ties indicate that this system is a good candidate of gap- 
less spin liquid. First, neutron scattering experiments re- 
vealed only short-range correlations of Ni spins even be- 
low ^20K. The correlation length saturates to £ ~ 20A, 
corresponding to seven times the inplane lattice con- 
stant 6 ) and very short. Second, magnetic specific heat 
shows a power-law dependence Cm oc T 2 in the tem- 



perature regime 0.35-4K, which signals gapless and lin- 
early dispersive modes of excitations. 6 ) Lastly, magnetic 
susceptibility approaches a finite value as temperature 
approaches OK, indicating the absence of a finite spin 
gap. 6 ) 

To clarify the origin of this gapless "spin liquid" be- 
havior in NiGa 2 S4, Tsunetsugu and Arikawa proposed a 
scenario of antiferro nematic order. This is equivalent to 
an antiferro spin quadrupolar (AFQ) order, where order 
parameters are quadrupole moments, Q^^ = ^(S^S^ + 

S»'S») - §S'(S' + 1 )<W- 8) The y investigated an S=l spin 
model with bilinear and biquadratic (BLBQ) couplings 
on the triangular lattice, defined as 

H=J2[JS i -S j +K(S i -S j ) 2 ]. (1.1) 

(id) 

First they showed this model has an AFQ order that fits 
the tripartite triangular lattice in the parameter region 
of < J < K using mean field approximation. This 
mean-field ground state is represented as 

|*MF> = J] \ S * = °)a,R ® \Sy = 0) B.R ® l& = 0) C ,R > 
R 

(1.2) 

where j labels three sublattices (A, B, and C) and 
\S a = 0) ■ R denotes the single-spin state with eigen- 
value of £ Q -operator (a — x,y, or z) at the j-sublattice 
site in the unit cell R. Then they studied low energy 
properties in the AFQ phase using a bosonic description 
of the excitation and obtained results qualitatively con- 
sistent with the three essential points in the experiments 
in NiGa 2 S4: (1) absence of magnetic long-range order; (2) 
nonvanishing susceptibility at zero temperature; (3) T 2 
behavior of the specific heat. Lauchli et al. independently 
studied the ferro quadrupolar (FQ) phase in the parame- 
ter region K < J < —2.5K of the same BLBQ model and 
obtained the results similar to the AFQ case. 9 ) Although 
these proposals are suggestive, we have to note that the 
origin of large effective biquadratic coupling, either posi- 
tive or negative, remains to be clarified. Further, a more 
direct identification of quadrupolar order is desired. 

Spin freezing is another unusual phenomenon observed 
in NiGa 2 S 4 . The magnetic susceptibility shows a kink at 
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Tj = 8.5K, and a small bifurcation between field cool- 
ing (FC) and zero field cooling (ZFC) values below T/. 6 ** 
Muon spin rotation (fi SR) experiments revealed qua- 
sistatic relaxation of Ni spin below Tf. 10 " 1 These results 
suggest a spin freezing transition at Tf. However, the 
characteristic of the spin frozen state below Tf is re- 
markably different from the case of canonical spin glass 
materials. Slow Ni-spin fluctuations with a time scale of 
/is exist and are rapidly suppressed upon application of 
magnetic field > lOmT. 10 - 1 In order to further investigate 
this spin-freezing transition, Nambu et al. studied the 
thermodynamic properties of Nii_ 2; Zn. l; Ga2S4, where Ni 
ions are partially replaced with nonmagnetic Zn ions. 11 ' 
They showed that the freezing temperature Tf decreases 
with increasing impurity concentration x. This is just op- 
posite to the case of canonical spin glass materials. It is 
also important that Tf scales with Weiss temperature, 
which is also the characteristic energy scale of the low 
temperature specific heat. 

The main purpose of this paper is to propose a novel 
mechanism of spin freezing that is consistent with the 
spin liquid behavior and the unconventional spin freezing 
in NiGa2S4. Assuming the existence of the AFQ order, we 
will introduce impurity magnetic moments in the system 
and study interaction between them mediated by low 
energy excitation in the AFQ order. We will then discuss 
a possibility of spin-freezing caused by this interaction. 

This paper is organized as follows. In §2, we will intro- 
duce a model for a single disorder, which induces mag- 
netic moments in otherwise AFQ ordered system. We 
also describe the strategy of our calculations. In §3 we 
will derive effective continuum models to describe low- 
energy excitation in the AFQ order. Using this, we will 
study the one-impurity problem in §4, to investigate the 
coupling between an individual magnetic impurity and 
low-energy excitations in the bulk. In §5, we will derive 
interactions between the impurity magnetic moments, 
mediated by the low-energy bulk excitation. Then we 
will discuss the possibility of spin-freezing caused by this 
interaction in §6. Finally §7 is a short summary. 

2. Model and Strategy 

We start with introducing a microscopic model of 
NiGa2S4 including magnetic disorder, which can explain 
observed spin freezing, in §2.1. We mainly study the case 
of T = and our basic assumption is that the system 
has the AFQ order. Then we describe our basic strategy 
of calculations in §2.2. 



Ordered spin _ _ _ 

quadripole M A ^ m W W £ W W 
momenta W^^^WW^W W 

Impurity 
mag i] etc 
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Fig. 1. Schematic diagram for the interaction between impurity 
magnetic moments. Dotted lines denote abstract interactions. 





Fig. 2. (a)Position of a sulfur site in NiGa2S4. (b)Triangle bond 
disorder, three sublattices, and lattice vectors ftj [i = 1,2,3) 
defined below eq. (4.15). 



2.1 Model 

One feature of the spin freezing in NiGa2S4 is that 
the freezing temperature Tf scales with the characteris- 
tic energy scale of low-temperature specific heat, which 
is ascribed to collective excitations in the AFQ order. 
This implies that the spin freezing has close relation- 
ship with the AFQ order. We consider that the freezing 
is induced by the interaction between disorder-induced 
magnetic moments, which is mediated by low energy ex- 
citations in the AFQ order, as schematically shown in 
Fig. 1. 

As a specific realization of this scenario, we construct 
a microscopic model as follows. First, for the pure bulk 
system to which impurities are to be introduced, we em- 
ploy the BLBQ model (1.1), which is a minimal model 
for the AFQ order. We investigate the parameter region 
< J < K, where the ground state has the AFQ order. 
Note that the bilinear interaction is rather simplified in 
this model. Large third-neighbor interaction is believed 
to exist in NiGa 2 S4, 6 ) but it is not included in our model. 
We use the model (1.1) for simplicity of the discussion. 
Second, we introduce disorders to the system. Here we 
focus on bond disorder, not site disorder, considering ex- 
perimental results. Small deficiency of sulfur concentra- 
tion drastically enhances the glassy behavior in hysteresis 
of the magnetic susceptibility. 12 ) Sulfur ions are located 
above or below the center of each triangular plaquette, 
as shown in Fig. 2(a), and their orbitals are in dominant 
exchange pathways. Disorder is induced in exchange in- 
teraction due to the vacancy of sulfur sites, and this is 
essential for the spin freezing behavior. Taking this into 
account, we introduce a different bilinear coupling for the 
bonds in a triangular plaquette at the impurity position. 
Thus, the model reads 

n.n. 

H = E [JS l -S J +K(S l -S J ) 2 ] 

n.n. 

+ J2 [J'Si-Sj+KiSi-Stf], (2.1) 

where D denotes randomly distributed triangular pla- 
quettes and an example configuration is shown in Fig. 



J. Phys. Soc. Jpn. 



Full Paper. 



Author Name 3 



2(b). Hereafter we call this individual triad of disorder 
bonds, simply impurity. We assume that the biquadratic 
coupling K in cq. (2.1) is not affected by impurities, since 
we focus on the behavior of the magnetic dipole moments 
induced by disordered exchange couplings, while the local 
variation in K does not yield significant results. Further, 
we study the ferromagnetic case J' < 0, in which the 
model is consistent with the scenario above, since three 
spins on an impurity plaquette tend to align and form 
an impurity magnetic moment as a whole. 

2.2 Basic Strategy 

Before starting calculations, it is useful to describe the 
framework and limitations of the present study. Our goal 
is to obtain the interaction between impurities in the 
AFQ ordered state. These interactions arise from the 
interference of the modulations of the AFQ order and 
they have two parts. The first contribution is related to 
the fact that each impurity deforms the nematic order 
pattern in the host locally around it, and it is given by 
the interference of this static order parameter deforma- 
tion between the impurities. The second contribution is 
mediated by the interactions of impurities and quantum 
excitations in the bulk. One impurity interacts with dif- 
ferent sets of excitations depending on impurity magnetic 
state. Those excitations propagate in the bulk and inter- 
act with another impurity, which is also dependent on 
the magnetic state of the second impurity, and this leads 
to impurity-impurity interactions. 

In the present study, we focus on the first contribution, 
i.e. the one given by static deformation of the nematic 
order due to impurity and neglect the contribution of dy- 
namical quantum excitations. This may be partially jus- 
tified by the fact that the AFQ order is stable in the pa- 
rameter region (0 < J < K) of the BLBQ model and the 
reduction of the nematic order parameter due to quan- 
tum fluctuations is quite small. 8 ' 9 ' This does not exclude 
the possibility that quantum fluctuations play some es- 
sential role, but this problem is beyond the scope of this 
study and should be examined in the future. However, 
based on a heuristic argument, we expect that dynamical 
quantum effects also lead to impurity-impurity interac- 
tions with similar nature, and we will discuss this briefly 
at the end of §6. 

To describe static deformations of the nematic or- 
der, we employ a site-dependent mean field approxi- 
mation. The phase space is restricted to the subspacc 
of site-factorized wave functions I^mf) = Yli^i wnere 
ipi denotes a one-spin wave function at site i. Local 
nematic and magnetic order parameters are given by 
= (Hh( S i S i+ S J S iMi) ~ lS(S + l)^ and 
m i = (V'il'S'ilV'i) respectively, and the energy of the cor- 
responding configuration is given by (f UF | H | $ MF ). In 
this way, the energy is a functional of these local fields 
{Q} and {m}, and will construct a "classical" Hamilto- 
nian describing this energy cost and its continuum limit. 
This is a classical model because only static deforma- 
tions are considered there. Low energy configurations 
within this approximation are accompanied with long- 
wavelength distortions of the order parameters, and this 
distortion is referred to as excitation in the following. 



Within this framework, we will evaluate the energy and 
configuration of the ground state with two impurities. 
The result shows peculiar nature of impurity-impurity 
interaction, and this interaction is a key of novel type of 
spin freezing, which can describe the peculiarity of the 
spin freezing phenomena observed in NiGa2S4. 

3. Continuum Theory of the Bulk 

The bulk part of the model (2.1) behaves as a medium 
of the interaction between impurities and only low en- 
ergy excitations in the AFQ order play a significant role, 
while detailed lattice structure is not important. It justi- 
fies replacing the bulk part with an effective field theory 
describing low energy excitation, and let us derive the 
effective model in this section. There arc two kinds of 
excitations. One is nonmagnetic excitation correspond- 
ing to deformation of the order of spin quadrupole mo- 
ments. The other is magnetic excitation, which induces 
magnetic dipole moments. We introduce field variables 
describing these excitations and derive effective models 
up to the second order in these fields. Up to this order, 
nonmagnetic and magnetic excitations are decoupled. 

It is convenient to introduce the following time-reversal 
invariant basis for each site: 13 ' 

|1) + |T) 



\x) 



V2 



|*> = -i|0>, (3.1) 



where |1) , |0) , and |1) denote the cigenstates of S z oper- 
ator with eigenvalues 1, 0, and —1 respectively. A general 
one-spin wave function is represented as 

\d) = d x \x) + d y \y) + d z \z) , d e C 3 (3.2) 

and we define two real vectors as the real and imaginary 
parts of d: 

c3 



U + IV, 



u, v e 



(3.3) 



These vectors satify the normalization condition|ix| 2 + 
\v\ 2 = 1, and can also satisfy the orthogonality u ■ v = 
and the condition \u\ > \v\ by choosing an appropri- 
ate phase factor. We choose a local phase satisfying 
these relations in the rest of this paper. Using this rep- 
resentation, the expectation value of the Hamiltonian 
(1.1) with regard to the site-decoupled wave function 
I^mf) = lX is written as 



MF 



MF/ 



^ {(4J - 2K)[( Ui ■ Uj){vi -vj)- (ui ■ Vj )(vi ■ Uj)] 

(ij) 

+K[(u l ■ Uj) 2 + (vi ■ Vjf + (m ■ v/f + (vi ■ Uj) 2 ]} . 

(3.4) 

On the basis of this expression, we will derive the effective 
model for nonmagnetic excitation in §3.1 and then will 
turn to the effective model for magnetic excitation in 
§3.2. 

3. 1 Effective Model for Nonmagnetic Excitation 

Let us derive the effective model for nonmagnetic ex- 
citation first. This describes the energy of configurations 
under the condition that the magnetic moment m = 
at any site. We will show the effective model is the 0(4) 
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nonlinear-<7 model. Since the magnetic moment is given 
bym= (S) = 2u x v, the condition m = corresponds 
to |u| = 1 and v = 0. In this case one-spin state is char- 
acterized by vector it, referred to as director. Note that 
this representation is double- valued. Two states with di- 
rectors ±u differ only by a phase factor and hence cor- 
respond to the identical physical state. The Hamiltonian 
(3.4) becomes 



U{r)x[z] = X [C(r)], 



(3.8) 



H 



(ij 



«j) S 



(3.5) 



In the ground state of the disorder-free bulk system, 
directors in each sublattice are spatially uniform and or- 
thogonal between different sublattices. Nonmagnetic ex- 
citations mean a long- wavelength distortion of this set of 
ordered directors. Hence, as a local order parameter, we 
can use an orthogonal triad of unit vectors, which cor- 
responds to the directors at mutually nearest neighbor 
sites belonging to the three sublattices. Strictly speak- 
ing, the orthogonality of the directors between nearest 
neighbor sites can be slightly violated, but we can as- 
cribe this deviation to spatial variation of the local order 
parameter. We represent this triad as a SO(3), or al- 
most equivalently a SU(2) rotation operation, and here 
we choose the latter for convenience. This representation 
has multivalueness since the directors are double-valued 
and also each SO(3) matrix has double representation in 
SU(2). Although such redundancy plays an essential role 
for topological excitation, it is not relevant to the effects 
of spin-wave like excitation, which is studied in this pa- 
per. (We will summarize the properties of topological ex- 
citation in this system in Appendix B, and the important 
point is that the non-Abelian fundamental group of the 
AFQ order parameter implies a nontrivial merging rule 
of topological excitations. ) Therefore, we hereafter use 
the spin-1/2 representation of the SU(2) group, i.e. SU(2) 
matrix. Due to the locality and the SU(2) invariance, we 
can expect the effective model will have a following form, 

H* = y J drTr [VU\r) ■ VU(r)] , (3.6) 

with SU(2) field variable U(r) and stiffness constant 
J v . After a standard paramctrization, this model is 
mapped to the 0(4) nonlinear-cr model. Let us now derive 
this effective continuum model from the lattice model 
(3.5) and also determine the value of J ff . We construct 
the correspondence between the triad and the matrix 
as follows. First we define a double- valued unit vector 
£(aj»), 77(^2), or C( x i) depending on sublattice as 

£(xi) = ±Ui i e A, rf(xi) = ±Ui i € B, 

C{x i ) = ±u l zeC. (3.7) 

Under appropriate choices of signs, all vectors vary in 
space slowly compared with the lattice constant in the 
low energy sector. Then we extend these vectors by in- 
terpolation to those defined for continuous spatial coor- 
dinate r and impose the orthogonality among the three 
fields. Lastly we define the SU(2) field variable U(r) by 

U{r)x[x] = *[£(r)], U(r)x[y] = x[v(r)], 



where x = (1,0,0), y = (0,1,0), z = (0, 0, 1), and X [N] 
denotes the two component spinor representation of the 
spin 1 /2 coherent state 14 ) corresponding to the unit vec- 
tor N: 



X[N] = 



;os(f)exp(-^ 
in(f)exp(+^ 



(3.9) 



using polar coordinates as N = (sin#cosc/>,sin#sin0 
, cost9). It is known that a SU(2) matrix is written with 
four real parameters as 

3 

U(r) =7r (r) + 5^iT J -7r J -(r) 

3 

7r M eR, ( M = 0,l,2,3), 5>2=1, (3.10) 

where Tj's(j = 1,2,3) are Pauli matrices. Hereafter we 
regard 7r M 's (fj, = 0, 1, 2, 3) as local order parameters. 

Next we rewrite the bond Hamiltonian (3.5) in terms of 
the field variables 7r M (r). We can derive straightforwardly 
the representation of the continuated director vectors, 
£(r), rj(r), and C( r )> sucn as 

£ = (^O + n l ~ n 2 ~ ^h^{^l^2 - 7T07T3), 



2(7ro7r 2 + 7Ti7r 3 )) 



(3.11) 



Using these relations, and rewriting ir(xi) and it{xj) as 
7r, and tv + Stt respectively, we have interaction 

= 4K (-TT 3 dTr - 7r 2 (57ri + 7ri<57r 2 + tt 5it 3 ) 2 

= H AB (r,r + 5r), (3.12) 

between the sites i and j when i E A, j e B. We similarly 
define Hgc an d Hqa- 

Summing over all the site pairs and taking the con- 
tinuum limit, we obtain the 0(4) nonlinear-cr model as 
the effective model for nonmagnetic excitations, as we 
expected in eq. (3.6). 



1 



3 



J dr(VTv) 2 



E«2 

ju=0 



1. 



(3.13) 



where = (8/\/3)K. As shown by Polyakov and Wieg- 
mann, this model is identical to the phenomenologically 
introduced effective model (3.6), 15 ^ and the value of the 
coupling constant J v is determined explicitly. This clas- 
sical model describes long-wavclcngth distortion of the 
AFQ order, which physically corresponds to the energy 
cost of instantaneous deformation. 16 ) Note that the or- 
der parameter space SU(2) reflects the character of anti- 
fcrro order. Due to the orthogonality of directors between 
different sublattices, there are locally three independent 
directions of distortion of the order, which correspond 
to rotations of ordered moments around three axes in 
the spin space. On the other hand, in the ferro quadratic 
order where directors point to the same direction regard- 
less of sublattice, with the rotation around the director 
the order remains unchanged and the number of inde- 
pendent directions of distortion is two. Therefore we can 
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easily conclude the effective model for the ferro quadratic 
order is 0(3) nonlinear-cr model. 

Finally we introduce further simplification for the ef- 
fective model. We describe the ground state of the pure 
bulk as 7r = 1, 7r a = 0, (a = 1, 2,3). Even in the pres- 
ence of an impurity, the bulk region is close to this ground 
state. Hence we assume |7r | -C no ~ 1, (a = 1, 2, 3). Ex- 
panding the effective model (3.13) in terms of the field 
variables n a (a = 1,2,3) and preserving only the lowest 
order terms of n a , we obtain the three-component mass- 
less Gaussian model: 



\m 1 



(C) 



1 (C) 
2 TO 2 



\m 2 



(A) 



1 (-4) 1 (B 



and tp 3 = \vn\ = 



■in.. 



IT = \ J, 



3 r 

J dr(Vir a ) 



(3.14) 



Note that contribution of 7r is 0(n^) and therefore 
neglected here. We will use this effective model to study 
nonmagnetic excitations afterwards. This simplification 
of the model is not appropriate for investigating the ex- 
citations where ir— field configuration strongly deviates 
from that in the ground state. For example, at finite tem- 
peratures, the long range order is destroyed as manifested 
by Mermin- Wagner theorem, 17 ) and the n— field varies 
in space among the whole parameter space even without 
disorders. To deal with such a situation, we have to go 
back to the original model (3.13). 

3.2 Effective Model for Magnetic Excitation 

Now we turn to derive an effective model for magnetic 
excitation, which will be identified as a three-component 
massive Gaussian model. In contrast to its nonmagnetic 
partner, this describes the energy of configuration with- 
out nonmagnetic excitation. This condition means that 
the directions of the directors u do not change from those 
in the bulk ground state, while \u\ can be smaller than 
the bulk value, \u\ = 1. Hence we can impose 

ml/x (i € A), Uil/y (t g B), Uil/z (i e C). (3.15) 

First let us introduce field variables describing mag- 
netic excitation. Basically we can choose three com- 
ponents of magnetic moments as the field variables in 
concern, but we have to pay attention to the follow- 
ing two points. First, condition (3.15) yields a restric- 
tion of magnetic moments. For example, in the A sublat- 
tice, the magnetic moment is restricted in the YZ plane, 
, since m — 2u x v. Second, cou- 



th 



-(0, 



in 



pling between magnetic moments is ferromagnetic or an- 
tiferromagnetic, depending on the parameter J/K. It is 
clarified by transforming the total Hamiltonian (3.4) to 
the form 



(Hi,) 



J - ^K ) rrii ■ nij 



+ K{{u t ■ Uj) 2 + (vi ■ vff + (ui ■ Vj) 2 + (vi ■ Uj) 2 }. 

(3.16) 

The interaction between and rrij is ferromagnetic 
(antiferromagnetic) when < J < K/2 (K/2 < J < 
K). Taking these points into account, we choose field 
variables as follows. In the ferromagnetic region < J < 
K/2, the three variables are defined as ipi = = 



2 m l 

1 (A) 

2 3 



, ^2 



1 (C) 
l m 2 



, which represent 



(here the subscript denotes a component, not a site), 
which correspond to "uniform" magnetizations, as the 
field variables varying gradually in space. In contrast, in 
the antiferromagnetic region K/2 < J < K, the field 

variables are -01 ~ ~ 
— \rn^\ and ip 3 
" staggered" magnetizations 

Next we express the Hamiltonian in terms of the field 
variables ip and derive an effective model. Using the field 
variables defined above, the Hamiltonian (3.4) is written 
as: 

n.n. 

H= Yl {H l3 +H jk +H kl ) (3.17) 

ieA,jeB,kec 

Hij = -147 - 27^ 3 ,^3,i + Ktylj + 03 2 J etc. (3.18) 

As in the previous discussion for the nonmagnetic part, 
we take the continuum limit, and the result is the three- 
component Gaussian model with mass terms: 

1 3 

= =J,„ V / dr [(V^ a ) 2 + (M 2 V^] , (3.19) 



a=l 



where the coupling constants are given by 



(3.20) 



where Z denotes the lattice constant. Hereafter we take 
lo = 1 for simplicity. The effective model (3.19), is mas- 
sive, reflecting short-range nature of magnetic correla- 
tions of the AFQ order. Note that the "mass" van- 
ishes as we approach the boundaries of the AFQ phase, 
J = 0, or K. This singularity implies that the system 
becomes unstable against magnetic excitation, and man- 
ifests the onset of a magnetically ordered state. Actu- 
ally, the mean-field ground state has ferromagnetic order 
for J < 0, and 120 degree antiferromagnetic order for 
J > K.V 

4. One Impurity Problem 

In this paper, we want to calculate interaction between 
two impurities in the AFQ order, and this interaction is 
mediated by coupling of each impurity and bulk excita- 
tions. Therefore, our next task is to study the problem of 
single impurity and calculate its coupling constant of the 
interaction with nonmagnetic and magnetic excitations 
in the bulk studied in the last section. The most impor- 
tant result is that induced nonmagnetic excitation field 
shows a power-law decay in space, much more extended 
than magnetic excitation field. 

Our strategy is as follows. For the bulk part we use the 
continuum theory developed in the last section, while we 
use the original spin variables for the impurity part. We 
then derive the bulk-impurity coupling by evaluating the 
energy on the bonds connecting these two parts. 

At this point we explain our nomenclature concerning 
the sites and bonds. We denote sites in the impurity as 
"core" sites, and sites on the border of the bulk which 
are connected with core sites by bonds as "shell" sites. 
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Fig. 3. Simplification of the one-impurity Hamiltonian. (a) Black 
(white) circles denote core (shell) sites. (b)The bulk is replaced 
by a continuum media and its inner boundary is approximated 
by a circle. 



They are shown as black and white circles respectively in 
Fig. 3(a). We call the other sites "bulk" sites. Further we 
call bonds connecting core sites and shell sites (two core 
sites) "core-shell" ("core-core") bonds, which correspond 
to red (black) bonds in Fig. 3(a). Finally we call the other 
bonds "bulk" bonds, which are shown as gray bonds in 
the figure. 

First, we represent states on shell and bulk sites using 
the two vector field variables 7r and if> defined in the last 
section, instead of original spin wavef unctions. For core 
sites, we keep using original spin wavefunctions. Next, 
we evaluate the energy of bulk bonds using the effective 
massless and massive Gaussian models derived in the last 
section, whereas energy of core-core and core-shell bonds 
are evaluated for the original lattice model (1.1). In doing 
this, spin wavefunction at a shell site is determined by 
the field variables at its position. Then, we approximate 
the boundary of the bulk part, which has a polygon shape 
in the original lattice model, by a circle of radius tq, to 
make it easier to evaluate the energy of the bulk part. 
Finally, we consider only the dominant components of 
spatial fluctuations of the fields. 

For simplicity, we take the limit J' — » — oo i.e. the mag- 
netic interaction between core sites is ferromagnetic one 
of infinite strength. In this limit, spins on three core sites 
are completely polarized and described by the identical 
magnetic moment vector of unit length m, \m\ = 1. We 
change m and clarify the anisotropy of the magnetic mo- 
ment. We will show later, using numerical calculations, 
that the cases of finite J' (J' < 0) are qualitatively sim- 
ilar to this limit. 

The total Hamiltonian is 



H — ^bulk 



rri> 
"bulk 



#imp, (4.1) 



where H£ ulk and H^ ulk are the effective Hamiltonians 
for the nonmagnetic and magnetic excitations (3.14) and 
(3.19), and -£/i mp denotes the bond Hamiltonians on core- 
shell bonds. -ff; mp represents impurity-bulk interaction 
and is a function of the magnetic moment m of core sites 
and the field variables tt for nonmagnetic excitation and 
ip for magnetic excitation on shell sites. 

First we will derive i?i mp perturbatively with regard 
to tt and tp up to the second order in §4.1. Then, we 
will consider only the first order term within H{ mp and 



examine the magnetic anisotropy. At this order we can 
decouple the total Hamiltonian (4.1) into nonmagnetic 
part H v and magnetic part . We will treat them sep- 
arately in §4.2 and §4.3, respectively. Next we will further 
include the second order term and refine the anisotropy. 
Since we can not decouple the Hamiltonian as before, we 
will treat the whole Hamiltonian of this order in §4.3. 
Finally we will approach this problem numerically and 
compare the results with analytical ones in §4.4. 

4-1 Perturbative Expansion of the Impurity- Bulk Inter- 
action 

Following the strategies mentioned above, firstly we 
expand the impurity-bulk interaction -ff; mp with regard 
to the two field variables, tt and ip. An individual bond 
Hamiltonian H- lmp j, connecting the core site i and the 
shell site j, is written in terms of the impurity magnetic 
moment m, and the field variables on j site: 7Tj and tpj, 
up to the second order of TTj and ipj . Here we explain for 
the case where the site j belongs to A sublattice. 

If there is no nonmagnetic excitation at a shell site, we 
can derive the correspondence between the two vectors 
(it and v) and the field of magnetic excitation tp, from 
the definition of ■0. As mentioned before, nonmagnetic 
excitation corresponds to rotation of local nematic order 
in the spin space. The directions of principal axes are 
given by the representation of the directors, eq. (3.11), 
from which we can identify the rotation matrix. Thus 
we represent u and v, in terms of both magnetic and 
nonmagnetic excitation: 

M = (l-2(7r2+7T 2 2 )- l -{^l-CT^l), 

2(7Ti7r 2 - 7T 3 ), 2(ir 2 + ttitt 3 )) 



(4.2) 



V = (27T 3 V>3 - 2^20-^2, ^3 - 271-10-02, -C02 ~ 27TlV> 3 ) , 

(4.3) 

where the third and higher order terms are omitted and 
cr = +1 (-1) for J < K/2 (J > K/2). As for the 
core site i with magnetic moment m, we can readily 
obtain the value of u and i>, from the relations among 



these three vectors: \u\ 



1, u ■ v 



0, and 



m = 2u x v. Combining these results with the bond 
Hamiltonian in eq. (3.4), we obtain the representation of 
Himp,j{rn,Tz{xj),il){x :j )) for j e A. 



I I TTTT 

^impj' — 11 imp J 



H7 

imp,.? 



imp,j 



H; 



(4.4) 
(4.5) 



H L P ,j =(2J - K) [-am 2 Mxj) + m 3 M*i)] ( 4 -6) 
#im PJ =2^ {to? [^(x^-nKxj)] - [ei ah TO Q 7r fc (a; 3 )] 2 } 



m 



K 

T 



impj 



(4.7) 

{m 2 [ip 2 (x 3 ) -Vi(sCj)] - [eiabma-ip^Xj)} 2 } 

(4.8) 

(<7+l)/2 I \ i I \ i 
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■Kl€labm a llj b (Xj)] , (4.9) 

where e^fc is Levi-Civita tensor, and summation with re- 
gard to repeatedly appearing indices is implicitly taken. 
This derivation for the other sublattices is straightfor- 
ward and we obtain the results by changing the indices 
of m, 7T, and ip in a cyclic way: 1—^2, 2 — >■ 3, 3 — > 
1, (j € B), 1 -> 3, 2 -> 1, 3 2, (j G C). Finally the 
total contribution of all the shell sites is 

HL P = E H U,j (I = 7T^,n 2 ,^,^). (4.10) 

jGshclI 

4-. 2 First Order Effect of the Interaction between the 
Impurity and Nonmagnetic Excitation 
At this point we focus on the coupling between the 
impurity magnetic moment and nonmagnetic excitation 
in the bulk region. In the impurity-bulk interaction -ffi mp , 
we examine here only the first order term with regard to 
nonmagnetic excitation 7r, given by eq. (4.5). The total 
Hamiltonian becomes 

Z a=l J 

(4.11) 

where J denotes integration over the bulk region B' with 
a circular void shown in Fig. 3, i.e. J' dr = dr d6. 

As mentioned before, we consider only the dominant 
component of the spatial fluctuation of ir. This means as 
follows: As is clear from H£ ulk (3.14), the ground state 
configuration of 7r satisfies the Laplace equation V 2 7r a = 
0, (a = 1,2,3), and therefore can be expanded in the 
polar coordinate system as 

oo 

**(r, 0) = £ {ci a) r- n cos[n(9 - 9$)]} , (4.12) 

n=l 

where we impose the boundary condition 7r a (a;) — > 
with | a; | — ¥ oo. The dipole component c^ a ' becomes dom- 
inant as | a; | —¥ oo, and we neglect the other components 
c>n \ (n > 1). Hereafter we call this simply "dipole ap- 
proximation" . In this approximation, the field variable ir 
is given by 

7r a (r,0) = ^|^, (4.13) 

where we introduced a dipole moment //J and its value 
will be determined afterward. Recall that ro is the radius 
of the circular void of the bulk region. The bulk part of 
the ground state energy becomes 

Ku lk = ^E(^) 2 - ( 4 - 14 ) 

a 

Let us study the ground state of the total Hamilto- 
nian (4.11), in the dipole approximation. To be specific, 
we consider the case that the impurity is located at po- 
sition E in Fig. 4. The other cases will be summarized 
in Appendix A. By straightforward calculation, we find 
that the energy of the lattice part Hamiltonian has 



a simple form: 

3 

EL P = -fK E ™m b m c K a -/x- (4.15) 

a,6,c— 1 

where re a 's are lattice vectors: k% = (1/2, — v3/2), k 2 = 
(1/2,^/3/2), k 3 = (-1,0), shown in Fig. 2(b). 

The total energy is the sum of the bulk part (4.14) and 
the impurity part (4.15): 

E=^J2(»:f-~Kj2 l ^rn b rn c k a .»Z (4.16) 

a abc 

The dipole moments are determined by minimizing this 
energy with respect to /j,™ and the result is 

b,c 

where /j,q — 48iS>o/(77rJ 7r ). The ground-state energy is 
E£{m) = -ttJ^I) 2 (mi + m i 2 + m\ - l) . (4.18) 

This expression reveals that the impurity magnetic 
moment has an anisotropy. There are four easy axes 
^(1, ±1, ±1). Note that these directions should be con- 
sidered relative to the principle axes of the local AFQ 
order, not to the spatial directions of the triangular lat- 
tice. Since the Hamiltonian (2.1) is SU(2) invariant, we 
can interpret this anisotropy as a result of SU(2) sym- 
metry breaking in the AFQ phase. 

The anisotropy energy (4.18) can also be represented 

as 

ES(M) = -2irJM 2 (M% y + M 2 Z + M 2 X ), (4.19) 

using t 2 g part of spin quadrupole moment of each impu- 
rity spin 

M aP = l -{S a S p + S S a ) (a,f3 = x,y,z, a£P), 

(4.20) 

since M a p = m a mfi/2 is satisfied for a fully polarized 
spin with 5 = 1. We can naturally understand this by 
noticing that the nonmagnetic excitation ir linearly cou- 
ples to the spin quadrupole moment M a p in the first 
order coupling (4.5) considered here. 




Fig. 4. Six distinct possibilities of the position of the impurity. 
The position "i?" is chosen as reference, while the other cases 
may be reproduced by applying a point group operation shown 
in the figure. 
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Fig. 5. Further simplification of the one-impurity Hamiltonian. 
Couplings between core (black) and shell (white) sites are re- 
placed by the coupling between the impurity magnetic moment 
and the boundary A of the bulk region. 



We expect the results (4.18) and (4.19) also explain the 
magnetic anisotropy for finite J' < 0, although the pref- 
actors may be renormalized. In this case the magnetic 
moment m and the spin quadrupole moment M should 
be interpreted as the average among three impurity sites. 
This expectation is consistent with the numerical result, 
which will be shown in §4.5, that the anisotropy energy 
for finite J' < as a function of J/K is qualitatively 
similar to those for J' = — oo. 

Before closing this subsection, we show for later use 
that it is possible to modify the single-impurity Hamilto- 
nian (4.11) so that the exact ground state coincides with 
the result of the dipole approximation (4.13) and (4.17). 
Detailed structure of lattice part is expected to be irrel- 
evant after the dipole approximation and we modify the 
Hamiltonian such that the single impurity magnetic mo- 
ment interacts with a field on the boundary A, a circle 
of radius r* , of the bulk region, as shown in Fig. 5(b). 
Under this simplification, the total Hamiltonian can be 
written as 



n —-"bulk 



TT7T 

imp 



IT 71 
"bulk 



1 3 

■5*5: 

a=l 



dx{Vir a f 



H L P = / dxf (m,ir(x)) , 
J A 



(4.21) 
(4.22) 

(4.23) 



where f denotes the integration inside the bulk region 
B" shown in Fig. 5(b). We require that the ground state 
of this Hamiltonian for impurity magnetic moment m is 
given by the result of dipole approximation (4.13) and 
(4.17), which we denote as 7r°(m). 

We can uniquely determine the Hamiltonian which sat- 
isfies these conditions: 

H " = \ J « E f dx t v (*■« - 7r "( m ))] 2 + E o (™)> 

(4.24) 

where E (m) denotes the ground state energy. Clearly 
7r° is the ground state of this Hamiltonian, and we can 
transform this Hamiltonian into the form of eq. (4.21)- 



(4.24) by integration by parts. Therefore, we regard eq. 
(4.24) as the simplified single-impurity Hamiltonian. As 
a result, the impurity part H-^ is derived as 



imp 
3 ,.„ 



"imp = n "bulk 



1 f" 

■ [V7T°(m)-2V7T ] + E%(m) 

3 „// 

= - J-kYI \ dx ^Kl{m) ■ Vtt + ££(m) 

a=l •* 

(4.25) 



4-3 First Order Effect of the Interaction between the 
Impurity and Magnetic Excitation 
Next we consider the effect of the coupling with mag- 
netic excitation in the bulk region. Including only the 
first order coupling of this kind, given by eq. (4.6), the 
Hamiltonian for magnetic excitation becomes 

H ^ =H hm + H imp , 



iJ =-J 

-"bulk — ~ J i 



3 
o=l 



dx [(VV-a) 2 + fc/V>«] , (4.26) 



Just like the discussion above for nonmagnetic excita- 
tion, we consider only the dominant component of the 
magnetic excitation. The ground state configuration of 
tf> satisfies the Helmholtz equation V 2 ip a — h$ ip a = 0, 
and can be expanded in the polar coordinate as 



l/>«(r,0) = E c^K m {k^r)cos[m{e- 



g(") \ 
7 0.m) 



(4.27) 



ro=Q 



under the boundary condition ^ a (x) — » with | as | — > 
oo. 18 ^ Here K m {x) denotes the modified Bessel function 
of the second kind. We naturally expect that the angular 
momentum components with large m are not domi- 
nant and neglect cf [m > 1). We call this "monopole- 
dipole approximation" . We will justify this approxima- 
tion by numerical calculations in §3.D. As a result, the 
ground state field configuration if) is 

i) a {r,6) = qtK {k^r) + K 1 {k i ,r)-r-nt, (4-28) 



where we introduced a charge qf and a dipole moment 
fif. Further, the bulk part of the energy becomes 



-^buik - J i> E 



where we defined 



£ Q (x) = ttxK (x)Ki(x), 



(4.29) 



(4.30) 



Now we study the ground state of the total Hamilto- 
nian (4.26), using the monopole-dipole approximation. 
We deal with the ferromagnetic (J < K/2) and antifer- 
romagnetic (J > K/2) regions separately, and first focus 
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on the former. Here we focus the impurity located at 
"E" in Fig. 4, and show the results for the other cases in 
Appendix A. In this case the energy of the lattice part 



Hamiltonian Hf becomes 



imp 



J^Y, [4{-Kb(*tyri)+#i(A^r 2 )}m g2? 



a=l 



+ 



(— Kiik^n) + — K^r^X m a k' a ■ (if 



(4.31) 



where we defined k[ = (-V3/2, -1/2), k' 2 = (\/3/2, 
— 1/2), K3 = (0,1), and the radial coordinates of shell 
sites T\ — a/21/3 and r 2 = 2\/3/3. The total energy is 
the sum of the bulk part (4.29) and the impurity part 
(4.31), and given by 

E = J^J2 [£o(Vo) {(it it^a) 2 - {qtfml} 

a 

+£ 1 (k^r ) J2{(nt- l4m a k' a f - 



where we defined 



% 



V3 



Mo = 



£o(k^r 
1 



(K (k^n) + K (k^r 2 )) 



(4.32) 



(4.33) 



1 4 

4£ 1 (fc, / ,r ) \n v r 2 



Therefore, the charge and the dipole moment are deter- 
mined by minimizing this energy as 

9a = 90 m a, Ma = Mo m <X> (4-35) 

and the ground state energy is 



Ei{m) = - J, 



(4.36) 

Note that this ground state energy does not depend on 
the direction of the impurity magnetic moment m. It 
means that magnetic excitation does not contribute to 
the anisotropy energy of the impurity magnetic moment, 
in the first order coupling discussed here. 

Next we turn to the antiferromagnetic region (J > 
K/2). The lattice part of the energy becomes 



*t P = -^E(^i(Vi) 

a=l ^ 



H Ki{k^r 2 ) m a k a ■ 



(4.37) 



Combining this with the bulk part (4.30) yields the total 
energy 

e =JipY, [^(Vo) (it) 2 +£i(V„) 



(Mo - Ptm a k a f 



where we defined /2q = \Z3fJ% . Therefore, the charge and 



~ib \ 2 

Mo rn a 



(4.38) 



the dipole moment are given by 
eft — — 
and the ground state energy is 



/1q m a K a , 



Mo 



(4.39) 



(4.40) 



Note that this energy does not depend on the direction 
of impurity magnetic moment m, just as in the ferro- 
magnetic region. It is notable that the first order effect 
of magnetic excitation in the bulk region does not yield 
the anisotropy of impurity magnetic moment in either 
ferromagnetic or antiferromagnetic region. We will see in 
the next subsection, however, magnetic excitation does 
contribute to the anisotropy energy through second order 
coupling. 

Just like in the previous subsection, we can mod- 
ify the single-impurity Hamiltonian (4.26) so that the 
monopolc-dipole approximation is exact: 

a=l 

+kl[ij a -yj° a (m)] 2 }+Et(m), (4.41) 

where denotes the field configuration of the ground 
state in monopole-dipole approximation. The impurity 
part -ff^p is derived as 



(4-34) Ht p =\j* E /" dx { ' [V^° - 2V^] 



+^VM-2Va]}+<(m) 



J, 



a=l 
3 

v-E 

a=l 



(4.42) 

Finally let us calculate the total induced magnetic mo- 
ments 5m = ^2ig COTe (S), and the total squared induced 
magnetic moments Sm^ q = ^2i^. cole {S} 2 ■ In the ferroma- 
gentic region, these are derived by using magnetic exci- 
tation ip as 

16tt 



5m 



3^3 



s < =371 



dr%l> — 
drif) 2 



■^ Tn ^ Kl{k ^ ro) (443) 



= ^=7rrg {(qt) 2 [Kl(k^r ) - K^r )] 
2 



k^r 



Koikip^K^k^ro) 



(4.44) 



These values diverge as 1/J for J — > 0. Note that 
this divergence should not be literally taken. For small 
J, magnetic excitation tp on the sites near the impu- 
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rity are not negligible and perturbative approach with 
regard to ip breaks down. Therefore we anticipate actual 
divergence of induced magnetic moment 5m is weaker 
than 1/J. In the antiferromagnetic region, the present 
approach does not predict the value of the total induced 
magnetic moment 5m, since the field variable ip cor- 
responds to the staggered magnetization. On the other 
hand, total squared moments m sq can still be represented 
by the first line of eq. (4.44) and we obtain 

4 ^rl(~4f [tfg(V-o) - tf?(**r ) 



5m, 



sq 



3\/3 



k^r 



(4.45) 



which diverges as 1/(1 — Kj J) for J — > K. Again, ac- 
tual divergence is expected to be weaker. We will numer- 
ically investigate these two representations of induced 
magnetic moment in §4.6. 

4-4 Second Order Effect of the Interaction between the 
Impurity and Bulk Excitation 
In the previous argument, we showed that the mag- 
netic anisotropy emerges from the first order coupling to 
nonmagnetic excitation in the bulk region, and that the 
anisotropy is corner-cubic one represented by eq. (4.18). 
We will confirm this, by means of numerical calculations 
in the next subsection, but we will also find different 
types of anisotropy appear when J/K is small. They 
originate from higher-order effects of the coupling be- 
tween impurity magnetic moment and bulk excitations. 
In this subsection we examine the effects of the second 
order coupling and investigate the magnetic anisotropy. 
Adding the three second-order terms (4.7)-(4.9), the to- 
tal Hamiltonian is now given by 



iJ (2) =i?tuik 



^bulk + ^imp 



imp 



imp 



H, 



imp 



imp 



(4.46) 



where the bulk parts H£ ulk and H^ ulk are given in eq. 
(4.11) and (4.26), respectively. As for the bulk region, we 
have employed the dipole approximation for the nonmag- 
netic excitation 7r, and the monopole-dipole approxima- 
tion for the magnetic excitation if>. Here we introduce two 
additional approximations, in order to simplify the prob- 
lem. First, we consider only the lowest-order nonzero an- 
gular momentum component for the magnetic excitation 
xjj. As is clear from the previous results (4.35) and (4.39), 
it corresponds to the monopole approximation for the 
ferromagnetic region (J < K/2) and dipole approxima- 
tion for the antiferromagnetic region (J > K/2). Numer- 
ical calculations in the next subsection shows that these 
lowest-order moments are more than several times larger 
than the higher order moments. Second, we assume that 
the dipole moments for the field variables point to the 
direction obtained in the first order calculation. Again, 
numerical calculations will verify this approximation. As 
a result of these two approximations, the field variables 
are represented as 



7r a (r) 



ro 



Mr) =qtK Q (k^r) 



F.M. (4.48) 



1 



Mr) =#i(^r)-/£ ■ k a ■ r A.F.M., (4.49) 

with scalar variables pF a ,q^, an d Ma- Note that it is 
straightforward, though not shown here, to extend the 
following argument to more general situation that we do 
not employ these approximations. 

Now let us calculate the ground state energy of the 
Hamiltonian (4.46) within the approximations above, 
and we first focus on the ferromagnetic region (J < K/2). 
For the configuration (4.47)- (4.49), the energy becomes 

E{2) = E 9i- M u g j -2 ]T fi -gi (4.50) 

where fj and gi, (I = n, tp) are three-component vectors 
/tt =Cnfio (m 2 m 3 ,m 3 m 1 ,m 1 m 2 ) T (4.51) 



(4.52) 
(4.53) 

(4.54) 



U = c v9o (m 1 ,m 2 ,m 3 ) T 

94, = {qf,qf,qf) , 

Mjj's are three-by-three matrices 

{M^ij = c n 5 l: j + -^2(1 - 5ij)mim,j, 

(M^) rij = c^Sij + -0^,2(1 - 5ij)mimj, 



1 3 

(M^)^ = (M^ij = ^ \tijk\m kl (4.55) 



fe=i 



with coefficients 



369 



: - ; -^~ r o> C V = J tl> £ o (k^ro) , c T 2 = - — Krl, 



c^2 =2K jifo (k^ri) 2 + K (fc^r 2 ) 2 j 



(4.56) 
(4.57) 



c w4 , = - 2(2 J - K)r 



(k^n) + -^K (k^r 2 ) 



;\ 1 a K a • r, 



(4.47) 



(4.58) 

The ground state energy of the Hamiltonian (4.50) is 
readily derived as 

E (m) = - ]T fj ■ (M„ - tjj^fj (4.59) 

Here, the self-energy part is given as 

Stttt =M^M^M^, = M^M^M^ (4.60) 

=M^M^M^, = M^M-jM^ (4.61) 

and represents the effects of coupling between nonmag- 
netic and magnetic excitations. (We do not show the ex- 
plicit form of the ground state energy, since it is quite 
lengthy.) Later we investigate the dependency of the 
ground state energy on the coupling constants J/K. 

Next, we turn to the antiferromagnetic region (J > 
K/2). The ground state energy is similarly calculated 



J. Phys. Soc. Jpn. 



Full Paper 



Author Name 11 



and the result is given by replacing the variables in eqs. 
(4.59)-(4.61) by those for the antiferromagnetic region as 



C^2 



K 



Cn/, =2(2 J - K)r 



C0 = J^£ 1 (k^ J r ). 



13 



(4.62) 

j-ifi(fc^r 2 ) 

(4.63) 
(4.64) 



and q$ in cq. (4.52) and qf (a = 1,2,3) in cq. (4.54) 
should be replaced with /2q and , respectively. 

We show later in Fig. 8 (b) the ground state energy 
(4.59) as a function of J/K when the impurity magnetic 
moment is fixed along the representative directions. We 
can see that the energy now depends on J/K, unlike the 
results with the first order coupling (4.18) and (4.36). 
This is an effect of the second-order couplings c v % , c^a 
and <V0. In the parameter region < J/K < 0.077, 
(100) direction is an easy direction, i.e. the anisotropy is 
"face-cubic". We will compare this result with numerical 
results in the next subsection. 

4-5 Numerical Calculation 

We have shown analytically the appearance of the 
magnetic anisotropy and the configuration of the ground 
state, based on the effective field theory for the low- 
energy excitations in the bulk region. In this subsection 
we numerically solve this one impurity problem, in order 
to compare with the approximate analytical results in 
the previous subsections. 

We calculated the ground-state energy and spin con- 
figuration, by numerical calculations with finite-size clus- 
ters up to L — 40, including 5292 sites. The geometry of 
the cluster is depicted in Fig. 6. L labels the layer of 
sites away from the impurity bond triad placed at the 
origin, and it characterizes the cluster size. To minimize 
finite-size effects, we set the spin wavefunctions on the 
outer-boundary sites (shown by black circles in Fig. 6) to 
the bulk values (1.2). With this boundary condition, we 
minimize the energy (3.4) by optimizing wavefunctions 
of each spin. As for the ferromagnetic coupling between 
core sites J' < 0, we examined both of finite and infi- 



nite case. For the infinite case (J' = — oo), we fixed the 
direction of impurity magnetic moment and examined 
the ground state for each value of the impurity magnetic 
moment, in order to compare with the analytical results. 

First let us show the results for finite J'. Fig. 7 shows 
the phase diagram and each phase is characterized by 
the easy directions of the average magnetic moment of 
the three core sites: m = | X^igcore m «- We find that 
there are five phases with different easy axes and label 
these phases as shown in Table I. The (111) phase has a 
corner-cubic anisotropy, and covers the largest region of 
the coupling J/K, when | J'\ is sufficiently large. On the 
other hand, the (100) phase has a face-cubic anisotropy 
and appears near the J = line. Moreover, between these 
two phases we find the (110) phase, which has an edge- 
cubic anisotropy, and the least symmetric (abO) phase. 
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Fig. 6. The geometry of the finite size sample used in the numer- 
ical calculation. The case of size L = 2 is shown. States on the 
sites denoted by black circles are fixed to the mean-field ground 
state (1.2). 



Fig. 7. Phase diagram of the ground state in J-J' parameter 
space, (a) < J/K < 1, (b) < J/K < 0.1. 



Table I. Labeling of five phases and easy directions of impurity 
magnetic moment in each phase. 



label 


easy directions 


(000) 


(0,0,0) 


(111) 


^(±1.±1,±1) 


(110) 


^(±1,±1,0), -^(±1,0, ±1), ^=(0,±1,±1) 


(abO) 


(±a,±6,0), (±b,±o,0), (±a,0,±6), 
(±6,0,±a), (0,±a,±fe), (0,±6,±a) (a 2 + b 2 = 1) 


(100) 


(±1,0,0), (0,±1,0), (0,0, ±1) 
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Fig. 8. (a)(b)Energies of the (111) and (110) states measured 
from the energy of the (100) state in the limit J' = — oo. 
(a):numerical result, (b) analytical result. (c)(d) The ground 
state phase diagram in the limit J' = — oo. (c):numerical result, 
(d): analytical result. 



Second we show the results in the limit J' = — oo. Fig. 
8(a) and (c) show the energy difference between three 
representative phases and the phase diagram, respec- 
tively. The corresponding analytical results are shown 
in Fig. 8(b) and (d). The four phases observed at finite 
| J' |, survive in this J' = — oo limit, as shown in Fig. 
8 (c). It implies that the finiteness of \J'\ is not essen- 
tial for the presence of these anisotropies of the impu- 
rity magnetic moment. Further, the appearance of the 
two phases (100) and (111), which cover large regions 
in the phase diagram, qualitatively agrees with the an- 
alytical results shown in Fig. 8 (d), whereas the other 
two phases appear only in numerical calculation. The 
dependency of the magnetic anisotropy energy, shown in 
Fig. 8(a), qualitatively agrees with the analytical results 
in Fig. 8(b), although the difference becomes significant 
near J/K = or 1. This discrepancy can possibly be as- 
cribed to the breakdown of the perturbative treatment of 
the magnetic excitation. Near the onset of magnetically 
ordered phases (J/K = or 1), the system becomes sen- 
sitive to the magnetic impurity and the amplitude of the 
field for magnetic excitation increases. 

Then we examine the ground state spin configuration 
(Si) in the limit J' — — oo. As a typical example, we fix 
the impurity magnetic moment as m = ^=(1, 1, 1). Bc- 




Fig. 9. Numerical result of spatial distribution of induced mag- 
netic moments in the limit J' = — oo and m = ^=(1,1,1). 
(a) ferromagnetic case (J/K = 0.05), (b) antiferromagnetic case 
(J/K = 0.95). Projections to XY and XZ planes are shown in 
the left and right panel, respectively. 



fore detailed argument, we show induced magnetic mo- 
ments around the impurity in the ferromagnetic region 
(J/K = 0.05) in Fig. 9 (a), and in the antiferromagnetic 
region (J/K = 0.95) in (b). In both regions, induced 
magnetic moments form a complex noncoplanar pattern. 
This noncoplanarity can be understood as the result of 
antiferro quadrupolar order. As mentioned in §3.2, mag- 
netic moments are perpendicular to the directors, which 
are orthogonal between different sublattices. Note that 
this orthogonality is violated around the impurity, due 
to the nonmagnetic excitation. 

In order to compare the results with analytical ones, 
we converted the spin configuration to the field vari- 
ables for nonmagnetic and magnetic excitation, by solv- 
ing equations (4.2) and (4.3). The results are shown in 
Fig. 10. Nonmagnetic excitation tti around the impurity 
is plotted for two-dimensional space in Fig. 10(a) and 
also along one direction in Fig. 10(b). We can see that it 
shows dipole-like angle dependence, changing sign once 
around the impurity site, and its radial dependence de- 
cays as 1/r. This result verifies the dipole approxima- 
tion employed in §4.2 at least on the qualitative level. 
The direction of the principal axis is k% = (1/2, — v3/2) 
and coincides with that in the analytical result (4.17). 
Next we turn to the magnetic excitation i/>, and show 
it in Fig. 10 (c) and (d) for the ferromagnetic region 
(J/K = 0.1) and Fig. 10 (e)(f) for the antiferromag- 
netic region (J/K = 0.9). Its angular part has isotropic 
(dipole-like) in the ferromagnetic (antiferromagnetic) re- 
gion, and the radial dependence decays exponentially in 
both cases. These results verify the monopole-dipole ap- 
proximation we used in §4.1. Again, for the antiferromag- 
netic region, the direction of the principal axis is ki and 
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Fig. 10. Numerical results of spatial distribution of the field vari- 
ables for J' = — oo and m = ^j(b 1> !)• (a):Nonmagnetic exci- 
tation 7Ti for J/K=0.5. (b):Log-log plot of 7Ti along the arrow in 
(a). (c):Magnetic excitation for J/K=0.1. (d):Scmilog plot of 
ipl along the arrow in (c). (e):tpi for J/K=0.9. (f):Semilog plot 
of ipi along the arrow in (e). (g):Decay length of ip±. Two insets 
show log- log plots of 5i/> = l/Aty near J/K = and 1, together 
with analytical prediction (3.20). 
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Fig. 11. (a): Induced magnetic moment Sm measured along the 
direction of impurity magnetic moment Sm = m ■ Sm and to- 
tal squared induced moment im, q , in the limit J' = — oo and 
m = ^=(1,1,1). (b): Log-log plot in the ferromagnetic region, 
(c): Log-log plot in the antiferromagnetic region. Notice the hor- 
izontal axis is 1 — J/K. Cluster size is carefully chosen so that 
the finite siza effect is negligible (less than 1%). 



cal and analytical values. This discrepancy is ascribed to 
the coarse graining of the lattice model, which is justified 
only when the characteristic length scale £,/, of the theory 
is much larger than the lattice constant. 

In order to further examine the validity of the simplifi- 
cation of the angular part for tt and tp, we calculated the 
amplitude of each partial wave in eqs. (4.12) and (4.27) 
up to quadrupole (m = 2) component. The results are 
listed in Table II. We can see that the dominant com- 
ponent (i.e. (fi*(J = 0.1,0.5,0.9), q^(J = 0.1), ^{J = 
0.9)) is more than several times larger than the other 
components. This justifies our approximation that only 
single dominant component is taken into account. 

Finally let us examine induced magnetic moments out- 
side the impurity Sm. Fig. 11 shows total induced mag- 



agrees with that in the analytical result (4.39). In Fig. 
10(g) we present the decay length £^ of magnetic excita- 
tion ip, which is derived by fitting the radial dependence 
of ip to the exponential form ip( r ) ^ exp(— r / '£^) . In the 
regions near the phase boundary with ferromagnetic (an- 
tiferromagnetic) phase, < J/K < 0.1 (0.9 < J/K < 1), 
£,/, agrees well with its analytical value £/, = 1/k^p. In the 
intermediate region 0.1 < J/K < 0.9 with £^ < 1, on 
the other hand, there exist discrepancy between numeri- 



Table II. Angular momentum components of the field variables 
for m = -^(1,1,1). q,fi, and d denote the absolute values of 
monopole, dipole, and quadrupole component respectively. 





J/K 


0.1 


0.5 


0.9 


7T 


if 


0.086 


0.138 


0.125 




0.013 


0.003 


0.029 


i> 




0.529 




0.000 




0.113 




0.925 




0.028 




0.207 
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netic moment 5m and total squared induced magnetic 
moment <5m 2 q , denned in §4.3, as functions of J/K for 
(111) phase. In the ferromagnetic region, 5m is paral- 
lel to the impurity magnetic moment and 5m and <5m 2 q 
diverge with J — > 0. This divergence is slightly weaker 
than 1/J, as mentioned in §4.3. In the antifcrromagnctic 
region, is antiparallel to the impurity magnetic moment, 
as anticipated from antiferromagnetic interaction. Con- 
trary to the ferromagnetic region, the absolute value of 
induced magnetic moment remains finite as J — > K. The 
squared moment m sq diverges with J/K — > 1 and the 
divergence is weaker than analytically derived behavior, 
(1- J/K)- 1 . 

5. Interaction between Impurities 

In the previous section, we showed that magnetic and 
nonmagnetic excitations are induced around an impurity 
magnetic moment. When there exist multiple impurities, 
interference of these excitations causes interaction be- 
tween impurities. In this section, we study this impurity- 
impurity interaction, using both analytical and numeri- 
cal methods. We will show that nonmagnetic excitation 
yields long-range interaction with spatial anisotropy. 

In order to derive interaction between impurities, let 
us consider a system with two impurities with their mag- 
netic moments fixed, and calculate its ground state en- 
ergy. As in the previous section, we describe the bulk 
part by the continuum effective model, as shown in Fig. 
5. Regarding the impurity part H lm p, we consider only 
the first order coupling of impurity magnetic moment 
and bulk excitations. As a result we can decompose the 
Hamiltonian in the magnetic and nonmagnetic parts as 
H = H w ({■*}) + H^({ip}), and study the ground state 
of each part separately, just like the discussion in the 
previous section. 

For each single impurity, we take only the most domi- 
nant angular component of excitations, tt and rp, in the 
bulk region. It corresponds to the dipole approximation 
for the nonmagnetic excitation 7r. As for the magnetic 
excitation tp, the dominant mode is a dipole component 
in the antiferromagnetic region (J/K > 2), and a scalar 
component in the ferromagnetic region (J < K/2). We 
use the minimal models of impurity (4.24) and (4.41), for 
which these approximations are exact. In these models, 
we take the limit r* — > 0, where r* denotes the radius 
of the void of the bulk region shown in Fig. 5(b), since 
finiteness of r* is irrelevant to the asymptotic behavior of 
the interaction energy as the separation tends to infinity. 

We will analytically derive the interaction mediated by 
nonmagnetic and magnetic excitations in §5.1 and §5.2, 
respectively, and then proceed to numerical analysis in 
§5.3. 

5.1 Interaction Mediated by Nonmagnetic Excitation 

Let us derive the interaction energy between two im- 
purities mediated by nonmagnetic excitation in the bulk. 
Using the simplified one-impurity Hamiltonian (4.21), 
(4.22) and (4.25), the total Hamiltonian with two im- 
purities is given by 

TTTT _ TTTX i TT^{ a ) i TT^i^) 

H --"bulk + W imp + W imp 



+ E Q (m^)+E (m^), (5.1) 

where mW and 7T°W, (i = a,/3) denote the impurity 
magnetic moments and the ground state configuration 
in the presence of the individual impurity, respectively. 
i = a,f3 denotes the impurity index. Note that tv° it- 
self depends on the magnetic moment at the impurity, 
7T°M = 7r°(x; mW). We can define the interaction energy 
H? ut as the difference between the ground state energy 
of the whole system, which is calculated for eq. (5.1), 
and the sum of the two ground state energies of the one- 
impurity Hamiltonians: 



Hint = E s-s. - 



J, jd ■ Vtt" ; 



(5.2) 



Using the expression of the one-impurity ground state 
(4.13) and (4.17), the interaction energy is expressed in 
terms of thrcc-by-three matrices 



TJ7T 



-nJ n k Ir 



(x 



(x) 



5abm { a x) ", (x = a,/3) 



ab 



cos(26» + ?f) 





cos(2fl 






cos 29 



(5.3) 
(5.4) 

.)• 

(5.5) 



where the separation between the two impurities is ex- 
pressed as rp — r a = r a [i = r Q( g(cos#, sin 6*). The in- 
teraction is biquadratic with regard to impurity mag- 
netic moments m, and its radial part has r -2 depen- 
dence, whereas its angular part has dipole-dipolc like 
anisotropy. Note that nonzero elements of T matrix are 
nothing but t2 g part of spin quadrupolc moment defined 
by eq. (4.20). Therefore the interaction (5.3) can also be 
interpreted as interaction between quadrupolc moments 
of impurity spins. 

5.2 Interaction Mediated by Magnetic Excitation 

Second let us derive the interaction between impuri- 
ties mediated by magnetic excitations, using simplified 
one-impurity Hamiltonian (4.41). The result is different 
between the ferromagnetic and antiferromagnetic regions 
{J%K/2). 

In the ferromagnetic region (J < K/2), the interaction 
energy is 

H? nt = -2^J^(qt fK a {k^r a , 3 )m^ ■ m^\ (5.6) 

where K (x) denotes the modified Bcssel function of the 
second kind. The long-range asymptotic form is given by 



2tt 3 



(5.7) 
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where we used the asymptotic form of the modified Bessel 
function of the second kind K v {z) as z — > oo: 

J 4^ 2 - 1 

r z H h-- 

2z 8z 



1 



(b) 



(5.8) 



The expression (5.7) shows that the interaction is ferro- 
magnetic and isotropic in both of the spin space and the 
real space, and exponentially decays with the character- 
istic length scale of magnetic excitation, 

-1/2 

-II . (5.9) 



1 



K 

\K-2J\ 



In the antiferromagnetic region (J > K/2) , the inter- 
action energy is 

= nJ^i) 2 m^-[K 2 {k^r aP )f + K a {k^r aP )l] 

(5-10) 

where the form factor / is defined in eq. (5.5) Its long- 
range asymptotic form is 



m 



;^*) 2 e-^mW-[/-l]m 



(0) 



(5.H) 

This result shows that, the bilinear coupling of the mag- 
netic moments and the exponential decay of the radial 
part are common with the ferromagnetic region, but 
the interaction has anisotropy. Its angular part contains 
dipolc-dipolc like anisotropic terms inadition to isotropic 
and ferromagnetic term. 

5.3 Numerical Calculation 

We have analytically derived the interaction between 
impurity magnetic moments, making use of the contin- 
uum field theory of the low energy excitation in the bulk 
region and simplified coupling between each impurity 
and the bulk. Now we numerically calculate interaction 
energy, in order to check the validity of these simplifica- 
tions. 

We calculated the ground state energy of a finite-size 
system with two impurities of L = 40 with the fixed 
boundary condition as explained in §4.2. Here L de- 
notes the number of layers around the origin and im- 
purities are located symmetrically around the origin. 
Then we determined the interaction energy by subtract- 
ing from the ground-state energy one-impurity contri- 
butions, which were also derived numerically. We took 
the limit J' — — oo and fixed the impurity magnetic mo- 
ments, for proper comparison with the analytical results 
in the previous subsection. In addition, we consider the 
cases (m^ a \m^) = (m A ,m B ), (mA,— ibb), where 
m A = ^(1,1,1), m B = ^(1,-1,-1). 

The analytical results show that the interaction con- 
sists of the part mediated by nonmagnetic excitation 
H* nt , and the part mediated by magnetic excitation Hf nt . 
One remarkable difference between the two parts is that 
the latter part changes sign, when either spin is reversed, 
while the former part does not. Using this property, we 
can separate the nonmagnetic and magnetic parts of the 
interaction as follows: 

H Lt =\ {H- m t{m A ,m B ) + H int (m A ,-m B )} (5.12) 
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Fig. 12. Distance dependence of (a) nonmagnetic interaction 
H7L t , and (c)(e) magnetic interaction (c) corresponds to 

J/K = 0.1 in the ferromagnetic region and (e) is for J/K = 0.9 in 
the antiferromagnetic region. Corresponding dominant angular- 
moment components are plotted in log-log scale in (b) and in 
semi- log scale in (d) and (f). 



H tt = 2 { H int(m A ,m B ) - H int (m A , -m B )} (5.13) 

Further, to clarify the angular dependence of each part, 
we performed partial-wave decomposition as H^f^ = 

co + Em=i c ™ cos ( m ( 6 ~ 6 'o m) )) ■ 

Let us present the results of these analyses. First, we 
show the nonmagnetic interaction H[ nt at J/K = 0.5 in 
Fig. 12 (a) and (b), where dependence on separation r a p 
and radial dependences of the three large components 
are shown. We can see that the dipole component C2(r) 
decays as r -2 , and this is the dominant component. It 
agrees with the analytical result (5.3). In addition, the 
principal axis of the angular dependence coincides with 
that of the analytical result (5.3). Next, we turn to the 
magnetic interaction Hf nt at J/K = 0.1 in the ferromag- 
netic region shown in Fig. 12 (c) and (d). The monopole 
component co(r) decays exponentially with separation r, 
and is dominant, which agrees with the analytical result 
(5.7). Lastly, the magnetic interaction Hf nt at J/K = 0.9 
in the antiferromagnetic region is shown in Fig. 12 (e) 
and (f). The dominant component is now dipole one 
C2(r), and the monopole component cn(r) is of the same 
order as dipole one but slightly smaller. This agree with 
the analytical result (5.11), and the ratio of the dipole 
and monopole components is Cq/c 2 = 0.34, at r = 10, 
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close to the analytical result, 0.41. Again, the principal 
axes agree with that of the analytical result. 

These numerical data show that the analytical results 
describe the impurity-impurity interactions correctly on 
qualitative level. 

6. Discussion on Spin Freezing 

Let us now apply the results in the previous sec- 
tion to discussion on unusual spin freezing observed in 
NiGa 2 S4. First we discuss the possibility that the long 
range impurity-impurity interaction, which is mediated 
by nonmagnetic excitation in the bulk region, causes 
freezing of the impurity magnetic moments. Based on 
this scenario of novel spin freezing, second we discuss 
the anomalies in experiments in NiGa2S4 concerning spin 
freezing, particularly persistent spin dynamics below spin 
freezing temperature Tf and also scaling behavior of Tf. 

Several features shown in our calculations are impor- 
tant to realize unusual spin freezing. We have shown that 
the anisotropy of impurity magnetic moment and short- 
range and long-range impurity-impurity interactions are 
caused by magnetic and nonmagnetic excitations in the 
bulk region. The anisotropic long range interaction me- 
diated by nonmagnetic excitation is particularly relevant 
to the spin freezing. Let us consider the case where impu- 
rities are randomly located and their density is small but 
finite. In this case the spatial anisotropy of the impurity- 
impurity biquadratic interaction yields randomness and 
frustration, which are the origins of freezing phenom- 
ena. Two other points are important. First, that the bi- 
quadratic part of the impurity-impurity interactions are 
long-ranged and decays as a power law, r~ 2 . Second, 
the nematic order in the bulk induces spin anisotropy 
of impurity magnetic moments and effectively spins have 
only discrete degrees of freedom. Gandolfi et al. studied 
Ising spins on a hyper cubic lattice with random inter- 
actions, and proved that it has a thermodynamical spin 
glass order at any temperature, if the interaction is suf- 
ficiently long-ranged. 19 ^ r~ 2 dependence on two dimen- 
sional lattices satisfies this condition. Spin discreteness 
and long-range interactions are essential factors stabi- 
lizing spin glass order in their theory, and our system 
also shares these two points. Therefore it is reasonable 
to expect that our system has a similar thermodynam- 
ically stable glass order, although the effect of different 
discrete spin symmetries is not clear. It is an important 
future issue to prove this expectation. Instead, here let 
us point out one expected feature of spin freezing. The 
impurity-impurity interaction is biquadratic with regard 
to impurity magnetic moments, and therefore is invari- 
ant under the local spin inversion of each impurity mag- 
netic moment. One could alternatively say that degrees 
of freedom involved in the interaction are headless vec- 
tors, which are derived from ordinary Heisenberg spins by 
identifying their two tips. Hence, this interaction should 
cause the freezing of the headless vectors, rather than 
original magnetic dipoles. To put it simply, in this freez- 
ing state, each impurity magnetic moment fluctuates so 
that the local expectation value of the dipole moment 
vanishes, although it tends to be parallel or antiparal- 
lel to the local easy axis which is spontaneously chosen. 



However, it is necessary to note that the direction of this 
easy axis is not statically fixed. Throughout this study 
the coordinate system in the spin space is defined relative 
to the ordered quadrupole moment in the bulk region as 
a reference frame. At finite temperatures this reference 
frame is not static. It fluctuates in the time scale of the 
autocorrelation time of spin quadrupole moment. There- 
fore the easy axis of the impurity magnetic moment also 
fluctuates with the same time scale. 

Next we discuss on internal magnetic field below Tf, 
on the basis of the freezing mechanism just proposed. 
As we mentioned in §1, fiSR experiments revealed the 
presence of randomly distributed internal magnetic field, 
which fluctuates with a time scale of fis, and this dy- 
namics of the internal field is suppressed under magnetic 
field Hl > lOmT. Our scenario can explain the origin 
of this slow relaxation: that is flip of the impurity mag- 
netic moments along their easy axes. Note that this flip- 
ping process without energy cost is a unique result of the 
biquadratic form of impurity-impurity interaction. This 
idea can also explain the suppression of the dynamics 
by magnetic field. The local degeneracy of the two direc- 
tions of impurity magnetic moment, which correspond to 
the two tips of the easy axis, is lifted by external mag- 
netic field due to Zeeman energy. Note that not only the 
impurity magnetic moment but also induced magnetic 
moments around the impurity, participate in this flip- 
ping process. Magnetic moments are induced around an 
impurity, as we revealed in §4. The radius of this region 
is the magnetic correlation length, which is divergingly 
large near the phase boundary, and also experimentally 
determined as about seven times the lattice constant, as 
mentioned in §1. Therefore, it is likely that relatively few 
impurities, for example 1% of the bulk spins, induce the 
quasi-static magnetic moments and their slow dynamics 
on most sites, like observed experimentally. It is also a 
future task to investigate this possibility on more quan- 
titative ground. 

Then we present a possible explanation for the scaling 
behavior of the freezing temperature Tf. As mentioned 
in §1, Tf scales with the characteristic energy scale of the 
low temperature specific heat upon controlling nonmag- 
netic impurity concentration. It implies that the freez- 
ing occurs simultaneously with some kind of transition 
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Fig. 13. An example of configuration of spin quadrupole moments 
around a quaternion vortex excitation in the AFQ order. On the 
dotted circle, we can not consistently define the coordinate axes, 
which are parallel to the director on the sublattices. 
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which is related to the bulk order. One candidate is the 
vortex unbinding transition. Topologically stable defects 
exist in the AFQ order, as is shortly discussed in Ap- 
pendix B, and therefore we expect the existence of a 
vortex binding-unbinding transition. This transition in 
the bulk affects impurity-impurity interaction. The AFQ 
order sets a reference frame of spin coordinate. For ex- 
ample, we defined x-axis as the director on A sublattice. 
In doing this, we chose one of the two opposite direc- 
tions, both of which correspond to the same director. 
Such a choice is arbitrary, but can be done consistently 
in space. This situation remains unchanged even at fi- 
nite temperatures if the system is in the vortex binding 
phase. In this phase vortices exist as bound pairs, and 
we can neglect their existence, except for the renormal- 
ization of the coupling constant, as far as the long range 
behavior of the system is concerned. The situation drasti- 
cally changes, on the other hand, in the vortex unbinding 
phase, where free vortices exist. As shown in Fig. 13, we 
can not define the coordinate system consistently around 
a vortex. Therefore the origin of the long range impurity- 
impurity interaction essentially breaks down in the un- 
binding phase. Although the present study does not show 
explicitly yet, it is natural that the interaction becomes 
short-ranged and the characteristic length scale is the or- 
der of the average distance between free vortices. Since it 
is known that two dimensional systems with short-range 
random interactions have no thermodynamical glass or- 
der, 20 ) we may say that the freezing temperature 7/ in 
our system coincides with the vortex unbinding transi- 
tion temperature. This is consistent with the observed 
scaling behavior of X/, since 7/ and the energy scale of 
the low-temperature specific heat are determined by the 
characteristic energy of the low-energy excitation in the 
AFQ order. We need to investigate the finite tempera- 
ture problem of our system in order to substantiate this 
idea. 

Now we briefly discuss the effect of the quantum fluc- 
tuation. One way to include this is spin-wave like ap- 
proach starting with the deformed nematic order. We can 
construct a bosonic Hamiltonian in this vacuum with ex- 
tended Holstein-Primakoff transformation, which is anal- 
ogous to those adopted for the AFQ ordered system with- 
out impurities. 8 ) As stated in §2.2, the quantum correc- 
tion to the interaction energy arises from some exchange 
process of the bosonic excitations. Even without quan- 
titative investigation, we can predict basic characteris- 
tics of this quantum corrections. Reflecting the spon- 
taneous AFQ order, there are gapless Goldstone modes 
with linear dispersion. 8 ) This implies the boson exchange 
process yields long-range (power-law decay) interactions. 
Furthermore, near the gapless point the excitations have 
nonmagnetic character, 8 ) which indicates the dominant 
long-range interaction is nonmagnetic one, while mag- 
netic bilinear interactions would be subdominant. There- 
fore we expect the basic characteristics of the impurity- 
impurity interaction and also the discussion on spin freez- 
ing will not be changed even if we include the quantum 
corrections to the interaction. 

Finally we shortly comment on other open issues. 
The models in the present study, both of the bilinear- 



biquadratic model for the bulk part and the triangular 
bond disorder for the impurity part, are basically phe- 
nomenological. Therefore we have to verify these models 
on a microscopic point of view. 

7. Summary 

We have studied the effect of triangular-shaped ferro- 
magnetic bond disorder in the 5=1 bilinear-biquadratic 
model on the triangular lattice, in the parameter region 
where the antiferro quadrupolar order is realized. We 
have shown that coupling between impurity magnetic 
moment and magnetic and nonmagnetic excitations in 
the bulk yields several kinds of anisotropy of impurity 
magnetic moment, depending on the coupling constant 
of the bulk. We have also demonstrated the existence of 
biquadratic, spatially anisotropic and long-range inter- 
action between impurity magnetic moments, and deter- 
mined their effective coupling constants. This interaction 
is mediated by nonmagnetic excitation in the bulk. Based 
on these, we have presented the possibility of glass-like 
freezing of impurity magnetic moments due to this inter- 
action. This scenario can explain the unusual spin freez- 
ing observed in NiGa2S4 with persistent spin dynamics 
and the scaling behavior of freezing temperature with the 
energy scale of the bulk. 
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Appendix A: Extension to General Configura- 
tion of Impurities 

In §4 and §5 we have considered the case that an the 
impurity occupies the position E in Fig. 4. Here we show 
results for the other positions. One impurity problem is 
explained in A.l and impurity interactions are given in 
A. 2 for general cases of pair configuration. 

A.l One Impurity Problem 

Here we extend the calculation of the ground state field 
configuration around a single impurity in §4.2 and §4.3, 
when the impurity occupies a general position in Fig. 4. 

As for the nonmagnetic part, the ground state configu- 
ration of nonmagnetic excitation 7r is given by eq. (4.13) 
and (4.17), for the configuration E in Fig. 4. For exam- 
ple, in the case of C3, we can obtain the ground state 
by |7r rotation in real space from that in the case of 
E. This corresponds to the change of the lattice vectors 
as Ki 4 it2, K2 — > K3, K3 -> Ki. The other cases are 
also similarly obtained. It is natural to relate the config- 
urations of the impurity, to the elements of the trigonal 
point group C3 V = {E, C3, C^ 1 , <7i, 02, 03}. That is why 
we introduced the notation for the configuration. 

Taking advantage of this relationship, the general ex- 
pression for the dipole moment (4.17) is given as: 

6,c 
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where h £ C$ v denotes the impurity position and 
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Here D] x is a three-dimensional representation of C3 V : 



D 



E = 




1 

1 

1 

1 

1 

1 

3 1 ( 

1 ( 

3 1 



(A-3) 

The magnetic part can be calculated similarly and we 
show only results. Ground state configuration of mag- 
netic field i/> is given by eq. (4.28), and monopole and 
dipole moments are given as 



it =1o m o,, 



t4 



Mo m a*^a i 



q 4 ' 



=0, 



jj,Q m a k h a , 



(J < K/2) 
(J > K/2), 



where we defined 



£.'h 
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£.'h 



(A4) 
(A-5) 

(A-6) 



Dh denotes the three dimensional representation of C3 V 
introduced by eq. (A-3), and 07, denotes A 2 representa- 
tion of C^-.cfe = &C3 = a c~ x = a "\ = CTcr 2 = a <?3 
I. 

Note that the ground state energy is independent of 
the impurity position for both nonmagnetic and mag- 
netic parts. 

A. 2 Interaction between Impurities 

Using the results of the previous subsection, we can 
calculate the interaction energy for general cases of im- 
purity pair configuration. The impurity a occupies the 
position h a , while the partner impurity j3 occupies the 
position hp in a unit cell far away from the impurity a. 
As for the nonmagnetic part, the interaction energy is 



1 a(i 



Tr (t™ f hahfj T^) (A-7) 



h a hp 



2ir 2ir 
D ai cos(2# + — ) + D a2 cos(20 - — ) 



+£» CT3 cos2^] D T hfj 



) ab 



(A-8) 



where matrices T^ x > , (x = a, (3) are defined in eq. 
(5.4). We can see that principal axes of dipole-dipolc like 
anisotropy depend on the position of impurity pair, but 
the interaction has the same form as before. 

As for the magnetic part in the ferromagnetic region, 
interaction energy (5.6) is independent of impurity pair 
configuration. In the antiferromagnetic region, the inter- 



Table B-l. Class multiplication table for the quaternion group, 
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action energy (5.10) is replaced by the general expression 

H? nt =nMilt) 2 m^ ■ [K 2 (k^r aP )f hahfi 

+K {k 4 ,r a p)l}mW. (A-9) 

Again, difference of impurity pair configuration is re- 
flected only in the change of principal axes of dipole- 
dipole like anisotropy. 

Appendix B: Topological Excitation in Anti- 
ferro Quadrupolar Order 

In §3 we mentioned the presence of topological exci- 
tations in the AFQ order, and it is closely related to 
the mechanism of spin freezing transition presented in 
§6. Here we briefly summarize the properties of these 
topological excitations. Detailed general arguments on 
topological excitations may be found, for example, in 
the review by Mermin. 21 ) In our context, topological ex- 
citations mean defects in static configuration of the or- 
der parameter which is not removable by any continuous 
transformation of configuration. 22 ) 

Since the system of our concern is two-dimensional, 
relevant topological excitations are point defects. They 
can not be removed by continuous transformation of local 
order parameters inside a contour enclosing the defect. 
The simplest example for this is a vortex in ferromag- 
netically ordered XY spins and this is characterized by 
an integer winding number. Only defects with the same 
winding number can be continuously deformed to each 
other. In general ordered media, a simple winding num- 
ber is not sufficient to label a defect, but one can use 
the homotopy theory. For order parameter space R, one 
defines the fundamental group tt\ (R) , which is generally 
not Abelian, and a topological defect can be labeled by 
one of its conjugacy classes. 

Order parameter space of the AFQ state is more com- 
plicated than simple cases like ferromagnetic order of 
XY spins. This is generally defined as the coset group 
G/H, where G is the symmetry group that keeps the 
original Hamiltonian invariant, while H is its subgroup 
that keeps the ordered state invariant. The bilinear- 
biquadratic model (1.1) has a complete spherical symme- 
try in spin space and the time reversal symmetry There- 
fore G is SO (3) x Z 2 - The AFQ state is invariant under 
180° rotations about three mutually perpendicular axes, 
which form the dihedral group D2 of order 4, and also 
the time reversal operation. Therefore the order parame- 
ter space is R = (SO(3) x Z 2 )/(D 2 x Z 2 ) = SO{3)/D 2 or 
cquivalently SU(2)/Q, where Q is the quaternion group. 
It is order 8 and non- Abelian, and has a two-dimensional 
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representation: 

Q = {1, -1, i°x, ~icr x ,ia y , ~i<r yi ia z , -icr z }, (B-l) 

where a a {a — x, y, z) are the Pauli matrices. 

It is known that for a continuous and simply connected 
group G, with a discrete subgroup H, the fundamental 
group of G/H is nx(G/H) = H. 21 ^ In our case, G is the 
rotation group SU(2) and H is the quaternion group Q, 
and therefore the fundamental group of our order param- 
eter space is m(R = SU(2)/Q) = Q. Note that SU(2) is 
a simply connected group but SO(3) is not. Topological 
defects are thus classified by its five conjugacy classes 
as: 21 ) 

Co = {1}, C = {-1}, 

C x = {±ia x }, C y = {±ia y }, C z = {±iv z }, (B-2) 

including Co corresponding to no topological defect. Note 
that three conjugacy classes contain multiple elements of 
Q. This is the result of non-Abelian nature of the group 
Q. 

The non-Abelian nature of the fundamental group Q 
implies nontrivial merging rules of topological defects. 
Consider two defects and what happens if we merge 
them. The point is non-uniqueness of the defect class 
for the merged defect. To define the defect class for the 
merged ones, one needs a closed path which encloses 
these two. The answer follows the class multiplication 
table for ir%(R) = Q shown in Table B-l. Note that some 
products contain more than one classes, reflecting non- 
commutative nature of the group Q. For that case, the 
answer depends on the configurations of other defects, if 
present. Typical cases are shown in Fig. B-l and the two 
paths both enclose the same two defects but they may 
yield different classes for the merged defect. 




Fig. B-l. Two contours Ti and T2 enclosing two defects going 
through the opposite sides of the third defect. Defects are shown 
by black circles. 
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